Versuch einer Objektdefinition über Minima

Vielleicht ist es möglich Objekte aus dem IR-10.8-µm-Kanal über die Minima und davon ausgehenden Grenzen zu definieren.

Um das zu untersuchen, brauchen wir zunächst ein paar Pakete.

In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import ndimage as ndi
from skimage.feature import peak_local_max

import cv2 

import datetime as dt

import l15_msevi.msevi as msv

import sys
sys.path.append("/vols/talos/home/stephan/utils/tracking")
import tracking_common as tc
/vols/talos/local/anaconda2-5.0.0/lib/python2.7/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  from ._conv import register_converters as _register_converters

Beispieldaten laden

Wir laden für den 28.07.2012, 12:00 UTC ein Beispiel.

In [2]:
t = dt.datetime(2012,7,28,12,0)

s = msv.MSevi(time=t,chan_list=["IR_108"])
s.lonlat()
s.rad2bt()

ir108 = s.bt['IR_108']
Region suggests use of hdf file

Nun suchen wir nach den lokalen Minima im Bild.

In [18]:
minima = peak_local_max(-ir108, min_distance=2,threshold_abs=-273)
In [9]:
print minima
[[597  88]
 [596  57]
 [596  45]
 ...
 [  2 450]
 [  2 398]
 [  2 255]]
In [19]:
h = 14
fig,ax = plt.subplots(1,1,figsize=(h*1.61803,h))
ax.imshow(ir108,vmin=210,vmax=300,cmap='gray_r')

x = [m[1] for m in minima]
y = [m[0] for m in minima]

ax.scatter(x,y,marker='+',color='blue')
Out[19]:
<matplotlib.collections.PathCollection at 0x7f9b7c3b1fd0>

Das sind ziemlich viele Minima. Es gibt mehrere Möglichkeiten deren Anzahl zu reduzieren:

  1. das Eingangsbild glätten und maskieren
  2. den Mindestabstand vergrößeren
  3. den Schwellwert strenger wählen.

Das probieren wir im Folgenden mal aus.

Glättung und Maskierung des Eingangsbildes

Für die Glättung gibt es mehrere Möglichkeiten. Hier probieren wir Folgende aus:

  1. Mittelwertglättung
  2. Gauß'scher Filter
  3. Medianglättung
  4. bilateraler Filter

All diese Verfahren sind abhängig von der Filterbreite. Was ist ein guter Wert dafür?

In [4]:
ir108_transformed = (tc.scale_array_min_max(ir108)*255).astype("uint8")

mean_blur = cv2.blur(ir108_transformed,(5,5))
gauss_blur = cv2.GaussianBlur(ir108_transformed,(5,5),0)
median_blur = cv2.medianBlur(ir108_transformed,5)
bilateral_blur = cv2.bilateralFilter(ir108_transformed,15,15,15) 
In [49]:
plt.imshow(ir108_transformed)
plt.colorbar()
Out[49]:
<matplotlib.colorbar.Colorbar at 0x7f9b7d4b1610>

Darstellung der Glättungsergebnisse

In [53]:
fig,ax = plt.subplots(2,2,figsize=(14,14))
ax[0,0].imshow(mean_blur,vmin=0,vmax=255,cmap='gray_r')
ax[0,0].set_title(u"Mittelwertglättung")
ax[0,1].imshow(gauss_blur,vmin=0,vmax=255,cmap='gray_r')
ax[0,1].set_title(u"Gaußglättung")
ax[1,0].imshow(median_blur,vmin=0,vmax=255,cmap='gray_r')
ax[1,0].set_title(u"Medianglättung")
ax[1,1].imshow(bilateral_blur,vmin=0,vmax=255,cmap='gray_r')
ax[1,1].set_title(u"bilateraler Filter")
Out[53]:
Text(0.5,1,u'bilateraler Filter')

Die Glättung mit dem bilateralen Filter sieht subjektiv am besten aus. Das Feld ist relativ stark geglättet, aber die Kanten sind noch gut erhalten.

In [115]:
thresh = (273 - np.min(ir108)) / (np.max(ir108) - np.min(ir108))*255
In [116]:
thresh
Out[116]:
140.29909227709862
In [117]:
minima = peak_local_max(-bilateral_blur, min_distance=2,threshold_abs=thresh)
In [118]:
h = 14
fig,ax = plt.subplots(1,1,figsize=(h*1.61803,h))
ax.imshow(ir108,vmin=210,vmax=300,cmap='gray_r')

x = [m[1] for m in minima]
y = [m[0] for m in minima]

ax.scatter(x,y,marker='+',color='blue')
Out[118]:
<matplotlib.collections.PathCollection at 0x7f9b7dc54310>

Das sieht schon besser aus. Nur die Häufungen von Minima sind seltsam.

Vergrößerung des Mindestabstandes

In [111]:
minima = peak_local_max(-ir108, min_distance=5,threshold_abs=-273)
In [112]:
h = 14
fig,ax = plt.subplots(1,1,figsize=(h*1.61803,h))
ax.imshow(ir108,vmin=210,vmax=300,cmap='gray_r')

x = [m[1] for m in minima]
y = [m[0] for m in minima]

ax.scatter(x,y,marker='+',color='blue')
Out[112]:
<matplotlib.collections.PathCollection at 0x7f9b7d9a3890>

Das sieht besser aus.

Strengere Schwellwerte

In [113]:
minima = peak_local_max(-ir108, min_distance=2,threshold_abs=-250)
In [114]:
h = 14
fig,ax = plt.subplots(1,1,figsize=(h*1.61803,h))
ax.imshow(ir108,vmin=210,vmax=300,cmap='gray_r')

x = [m[1] for m in minima]
y = [m[0] for m in minima]

ax.scatter(x,y,marker='+',color='blue')
Out[114]:
<matplotlib.collections.PathCollection at 0x7f9b7bdf9910>

Das sind immernoch viele Minima, aber sinnvoller verteilt.

Kombination aus allem

In [5]:
th = 250

thresh = (th - np.min(ir108)) / (np.max(ir108) - np.min(ir108))*255

test = 255 - bilateral_blur
thr = 255 - thresh

minima = peak_local_max(test, min_distance=15,threshold_abs=thr)

#from skimage.morphology import extrema
#minima = extrema.h_minima(bilateral_blur,30)
In [6]:
minima
Out[6]:
array([[471, 317],
       [471, 316],
       [404,  43],
       [404,  42],
       [404,  41],
       [404,  40],
       [403,  42],
       [403,  41],
       [398,  61],
       [398,  60],
       [392, 119],
       [392, 118],
       [391, 231],
       [391, 230],
       [390, 231],
       [389, 124],
       [378, 151],
       [370, 171],
       [370, 170],
       [370, 169],
       [369, 170],
       [287, 446],
       [280, 265],
       [280, 264],
       [279, 287],
       [278, 288],
       [278, 287],
       [277, 287],
       [276, 399],
       [276, 353],
       [276, 352],
       [275, 287],
       [274, 288],
       [274, 287],
       [273, 288],
       [270, 422],
       [269, 422],
       [269, 421],
       [259, 385],
       [259,  25],
       [259,  24],
       [258,  27],
       [254,  60],
       [251, 474],
       [251, 473],
       [245, 570],
       [245, 569],
       [245, 568],
       [245, 567],
       [244, 570],
       [243, 438],
       [243, 437],
       [241,  99],
       [240, 102],
       [240, 101],
       [235, 142],
       [230, 543],
       [230, 542],
       [230, 541],
       [230, 540],
       [230, 460],
       [229, 544],
       [229, 543],
       [229, 542],
       [229, 541],
       [229, 461],
       [229, 460],
       [229, 160],
       [229, 159],
       [229, 158],
       [228, 544],
       [222, 373],
       [214, 493],
       [213, 493],
       [209, 168],
       [208, 169],
       [208, 168],
       [208, 167],
       [207, 172],
       [207, 171],
       [207, 170],
       [207, 169],
       [203, 518],
       [203, 183],
       [200, 573],
       [200, 572],
       [199, 573],
       [199, 572],
       [196, 612],
       [178, 701],
       [172, 545],
       [172, 544],
       [171, 545],
       [171, 544],
       [170, 427],
       [170, 426],
       [169, 489],
       [169, 488],
       [169, 487],
       [169, 430],
       [169, 427],
       [169, 426],
       [168, 488],
       [168, 430],
       [167, 586],
       [156, 453],
       [153, 610],
       [153, 609],
       [153, 608],
       [152, 610],
       [152, 609],
       [152, 608],
       [152, 607],
       [151, 609],
       [151, 608],
       [151, 607],
       [151, 606],
       [151, 605],
       [151, 604],
       [147, 562],
       [147, 561],
       [133, 596],
       [120, 127],
       [119, 127],
       [117, 705],
       [116, 705],
       [104, 166],
       [100, 129],
       [ 99, 130],
       [ 98, 250],
       [ 96, 133],
       [ 96, 132],
       [ 92, 525],
       [ 92, 524],
       [ 92, 523],
       [ 92, 522],
       [ 91, 618],
       [ 91, 526],
       [ 91, 525],
       [ 91, 522],
       [ 85, 587],
       [ 85, 586],
       [ 82, 573],
       [ 82, 572],
       [ 80, 482],
       [ 80, 481],
       [ 79, 482],
       [ 79, 481],
       [ 79, 480],
       [ 78, 480],
       [ 78, 479],
       [ 77, 480],
       [ 77, 479],
       [ 77, 478],
       [ 76, 479],
       [ 76, 478],
       [ 76, 477],
       [ 73, 609],
       [ 72, 272],
       [ 69, 165],
       [ 69, 164],
       [ 68, 275],
       [ 68, 168],
       [ 68, 166],
       [ 68, 165],
       [ 68, 164],
       [ 68, 163],
       [ 67, 166],
       [ 67, 165],
       [ 63, 523],
       [ 63, 522],
       [ 63, 521],
       [ 63, 520],
       [ 62, 522],
       [ 62, 521],
       [ 62, 520],
       [ 61, 522],
       [ 61, 521],
       [ 61, 520],
       [ 60, 521],
       [ 60, 520],
       [ 60, 220],
       [ 60, 219],
       [ 58, 225],
       [ 58, 224],
       [ 57, 225],
       [ 55, 566],
       [ 54, 710],
       [ 54, 709],
       [ 53, 592],
       [ 52, 295],
       [ 49, 483],
       [ 49, 482],
       [ 48, 484],
       [ 44, 621],
       [ 43, 620],
       [ 43, 619],
       [ 42, 619],
       [ 42, 618],
       [ 31, 251],
       [ 26, 555],
       [ 25, 558],
       [ 21, 521],
       [ 21, 520],
       [ 20, 493],
       [ 20, 492],
       [ 20, 489],
       [ 20, 437],
       [ 19, 493],
       [ 19, 492],
       [ 19, 490],
       [ 19, 489],
       [ 18, 493],
       [ 18, 492],
       [ 18, 491],
       [ 18, 490],
       [ 18, 215],
       [ 18, 214],
       [ 17, 217],
       [ 17, 216],
       [ 16, 437],
       [ 15, 437]])
In [7]:
test = 255 - bilateral_blur
In [12]:
from mpl_toolkits.axes_grid1 import make_axes_locatable

def colourbar(mappable):
    ax = mappable.axes
    fig = ax.figure
    divider = make_axes_locatable(ax)
    cax = divider.append_axes("right", size="5%", pad=0.05)
    return fig.colorbar(mappable, cax=cax)
In [13]:
h = 14
fig,ax = plt.subplots(1,1,figsize=(h*1.61803,h))
#ax.imshow(np.ma.masked_less(minima,1))
test_plot = ax.imshow(test,vmin=0,vmax=255,cmap='gray')
colourbar(test_plot)
Out[13]:
<matplotlib.colorbar.Colorbar at 0x7f2f4dacbc50>
In [9]:
print(np.max(ir108) - np.min(ir108))
print(th - np.min(ir108))  
print thresh
109.76179048091547
37.39011596541775
86.86519716384645
In [14]:
h = 14
fig,ax = plt.subplots(1,1,figsize=(h*1.61803,h))
ax.imshow(ir108,vmin=210,vmax=300,cmap='gray_r')

x = [m[1] for m in minima]
y = [m[0] for m in minima]

ax.scatter(x,y,marker='+',color='blue')
Out[14]:
<matplotlib.collections.PathCollection at 0x7f2f40185d10>
In [15]:
h = 14
fig,ax = plt.subplots(1,1,figsize=(h*1.61803,h))
test_plot = ax.imshow(np.ma.masked_less(test,thr))
colourbar(test_plot)
Out[15]:
<matplotlib.colorbar.Colorbar at 0x7f2f401c1210>
In [190]:
h = 14
fig,ax = plt.subplots(1,1,figsize=(h*1.61803,h))
test_plot = ax.imshow(np.ma.masked_greater(ir108,250))
colourbar(test_plot)
Out[190]:
<matplotlib.colorbar.Colorbar at 0x7f9b37bc5850>

Das sieht schon recht gut aus. Bloß dies Häufung von Minima ist seltsam. Das kann man wahrschenlich folgendermaßen beseitigen:

  1. Matrix mit Nullen erzeugen
  2. Punkte die ein Minimum haben auf 1 setzen
  3. diese Matrix markieren
  4. Massenschwerpunkte der verschiedenen Objekte bestimmen und als neue Minima betrachten
In [16]:
from scipy import ndimage as ndi

new_minima = np.zeros_like(ir108)

for m in minima:
    new_minima[m[0],m[1]] = 1
    
lab,nlab = ndi.measurements.label(new_minima)

new_mins = ndi.center_of_mass(new_minima,lab,np.arange(1,nlab))
In [17]:
new_mins
Out[17]:
[(15.5, 437.0),
 (17.0, 216.5),
 (18.0, 214.5),
 (18.90909090909091, 491.27272727272725),
 (20.0, 437.0),
 (21.0, 520.5),
 (25.0, 558.0),
 (26.0, 555.0),
 (31.0, 251.0),
 (42.5, 619.0),
 (44.0, 621.0),
 (48.0, 484.0),
 (49.0, 482.5),
 (52.0, 295.0),
 (53.0, 592.0),
 (54.0, 709.5),
 (55.0, 566.0),
 (57.666666666666664, 224.66666666666666),
 (60.0, 219.5),
 (61.75, 521.0833333333334),
 (68.0, 164.75),
 (68.0, 168.0),
 (68.0, 275.0),
 (72.0, 272.0),
 (73.0, 609.0),
 (77.84615384615384, 479.6923076923077),
 (82.0, 572.5),
 (85.0, 586.5),
 (91.57142857142857, 523.8571428571429),
 (91.0, 618.0),
 (96.0, 132.5),
 (98.0, 250.0),
 (99.0, 130.0),
 (100.0, 129.0),
 (104.0, 166.0),
 (116.5, 705.0),
 (119.5, 127.0),
 (133.0, 596.0),
 (147.0, 561.5),
 (151.76923076923077, 607.6923076923077),
 (156.0, 453.0),
 (167.0, 586.0),
 (168.5, 430.0),
 (168.75, 488.0),
 (169.5, 426.5),
 (171.5, 544.5),
 (178.0, 701.0),
 (196.0, 612.0),
 (199.5, 572.5),
 (203.0, 183.0),
 (203.0, 518.0),
 (207.625, 169.25),
 (213.5, 493.0),
 (222.0, 373.0),
 (229.33333333333334, 542.2222222222222),
 (229.0, 159.0),
 (229.33333333333334, 460.3333333333333),
 (235.0, 142.0),
 (240.0, 101.5),
 (241.0, 99.0),
 (243.0, 437.5),
 (244.8, 568.8),
 (251.0, 473.5),
 (254.0, 60.0),
 (258.0, 27.0),
 (259.0, 24.5),
 (259.0, 385.0),
 (269.3333333333333, 421.6666666666667),
 (274.0, 287.5),
 (276.0, 352.5),
 (276.0, 399.0),
 (278.0, 287.25),
 (280.0, 264.5),
 (287.0, 446.0),
 (369.75, 170.0),
 (378.0, 151.0),
 (389.0, 124.0),
 (390.6666666666667, 230.66666666666666),
 (392.0, 118.5),
 (398.0, 60.5),
 (403.6666666666667, 41.5)]
In [19]:
h = 14
fig,ax = plt.subplots(1,1,figsize=(h*1.61803,h))
ir_plot = ax.imshow(ir108,vmin=210,vmax=300,cmap='gray_r')
colourbar(ir_plot)

x = [m[1] for m in new_mins]
y = [m[0] for m in new_mins]

ax.scatter(x,y,marker='+',color='blue')
Out[19]:
<matplotlib.collections.PathCollection at 0x7f2f4daef110>

Das sieht als repräsentativer für die Objekte aund Startpunkte für die Objektdefinition ziemlich gut aus. Es sind Minima von nichtkonvektiven Wolken dabei, aber um diese zu entfernen, ist wahrscheinlich eine Wolkentypisierung nötig.

In [27]:
ct = msv.MSGtools.get_nwcsaf_prod('CT',t,calibrate=True)
In [28]:
ct = np.ma.masked_where(np.logical_or(ct<5, ct>10),ct)
In [29]:
h = 14
fig,ax = plt.subplots(1,1,figsize=(h*1.61803,h))
ct_plot = ax.imshow(ct)
colourbar(ct_plot)
Out[29]:
<matplotlib.colorbar.Colorbar at 0x7f2f4d909090>